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Abstract 

We propose a 2-d computational model-system comprising a mixture of spheres 
and the objects of some other shapes, interacting via the Lennard-Jones potential. 
We propose a reliable and efficient numerical algorithm to obtain void statistics. 
The void distribution, in turn, determines the selective permeability across the 
system and bears a remarkable similarity with features reported in certain biological 
experiments. 
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1 Introduction 



The problem of packing of spheres plays a major role in the modeling of many 
physical systems and has been studied for more than four decades. Some of the 
early examples [1-3] of the computer simulations of hard sphere liquids sug- 
gest the existence of a first order freezing transition. The problem of packing 
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of spheres in two and three dimensions is of great interest. Recent investiga- 
tions of such systems have focused on the study of the statistical geometry 
of the dense sphere packing. Such studies are important in the understand- 
ing of physical properties of many systems, composed of a large number of 
particles [4-13]. 

In this context we pose a question, with the motivation of studying the trans- 
port across a two dimensional structure of packed circular disks (membrane) , 
how does the packing change when the membrane is doped with objects of 
various shapes and sizes (e.g. spheres arranged rigidly in the form of rods of 
different lengths, L, T, X shapes etc. See Fig. 1) ? In particular we investigate 
the effect of these shapes on the distribution of "voids". The "anisotropy" in 
the interaction potential appears to play a key role in the induction of large 
voids. 

As pointed out by Sastri et. al. [11], no algorithm is available to compute 
void statistics for the packing of shapes other than spheres. In this paper we 
propose a simple numerical algorithm to compute void statistics. Unlike a 
probabilistic algorithm (Monte Carlo), our algorithm is based on digitization 
and cell counting. 

The paper is organized as follows. In Sec. 2, we describe the model system. 
A definition of "void" and an algorithm to compute void statistics is given 
in Sec. 3. The results of numerical simulations and their relevance in lipid 
biomembranes is discussed in Sec. 4 We summarize the paper in Sec. 5. 



2 The model system 

The configuration space of the model system (membrane) is considered as 
a two dimensional space with periodic (toroidal) boundary conditions. The 
constituents of the membrane are disks and dopants. 



2.1 The basic model 

We consider a membrane made up of only circular disks interacting pairwise 
via the Lennard- Jones potential: 
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where, is the distance between the centers of the i and j disks, a deter- 
mines the range of hard core part in the potential and e signifies the depth of 
the attractive part. We choose the number of disks such that the area occu- 
pied by these disks is around 70%, which is less than that of the close-packed 
structure but still large enough to produce some closed voids. 

2.2 The model with impurities 

Further, we consider different shape anisotropic combinations (dopants) con- 
sisting of k number of circular disks. We treat each of these combinations as 
a single rigid cluster. Several such dopants (impurities) are considered. Fig. 
1 shows some of these impurities. The interaction between impurities and 
disks or other impurities is obtained by superposing the Lennard- Jones po- 
tential corresponding to each of the constituent disk in impurity. We consider 
a membrane with circular disks and impurities amounting to 10% of the total 
number of circular disks, such that the area occupied is still 70%. 

These membranes are brought to an equilibrium configuration by the Monte 
Carlo method [14] at a fixed temperature. Fig. 2 and Fig. 3 show typical equi- 
librium configurations of membrane without and with impurities respectively 
(The impurity in Fig. 3 is a rod shaped structure made up of five disks (Rod 5 ), 
in general Rod K for rod made up of k number of disks). In the simulation the 
temperature is so chosen that ksT < 4e, where ks is the Boltzmann constant. 
The equilibrium is confirmed by simulated annealing. 



3 Voids and an algorithm for void statistics 

Now, we introduce the notion of an "r-void" in a membrane which is suitable 
for the description of transport across membrane and further, propose an 
algorithm to compute statistical quantities such as the number of voids in the 
membrane, the void size distribution etc. 

We define an r-void as a closed area in a membrane devoid of disks or impu- 
rities, and big enough to accommodate a circular disk of radius r. Of course 
an r-void is also an r'-void if r' < r. 

3. 1 The algorithm to compute void statistics 

To compute the void statistics for r-voids, we increase the radii of the disks 
forming the membrane (including the disks in the impurities, without altering 
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the positions of the centers) by an amount r (See Fig. 4). Then we digitize the 
entire membrane on a suitably chosen grid. The choice of grid size depends 
on the required accuracy and the typical sizes of the voids. The digitization of 
circular disks is carried out by the Bressenham circle drawing algorithm [15], 
modified to incorporate periodic boundary conditions. The number of voids 
in the membrane are computed by flood filling [15] every closed void with a 
different color and then counting the number of colors. The sizes of various 
voids can be obtained by counting the number of grid-cells filled by the cor- 
responding color. The termination of flood fill algorithm is ensured since the 
voids are closed. In our case this condition is automatically fulfilled in view of 
periodic boundary conditions. 

The geometric algorithms involving Vorenoi polygons [11-13] are mathemat- 
ically satisfying and are expected to be accurate but would take much more 
computation time. On the other hand, as pointed in [11], the probabilistic 
algorithm is time efficient but requires a very large sample size while dealing 
with small voids. 

Our algorithm is quite efficient as well as suitable even when there are small 
voids in the membrane. We further note that the algorithm can be easily gen- 
eralized to higher dimensions. We expect that the efficiency of this algorithm 
can be further enhanced by the use of a multi-resolution adaptive grid. 



4 Results and Discussions 

The simulations were carried out for membranes of different compositions. Fig. 

5 shows the graphs of the number of r-voids as a function of r measured in 
units of the radius of the constituent disks. Curve (a) shows void distribution 
in absence of impurities. Curve (b) represents the void distribution in a mem- 
brane with rod shaped impurities made up of two disks (Rod2). Curves (c) and 
(d) show the void distribution with L shaped impurities made up of four disks 
(L 4 ) and rod like impurities made up of four disks (Rod 4 ) respectively. It is 
clear from the graph that the number of large voids increases with an increase 
in the anisotropy of the impurity. Even though L 4 and Rod 4 occupy the same 
area, Rod 4 being more anisotropic induces a larger number of big voids than 
L 4 . This fact can be clearly seen in Fig. 5, curves (c) and (d). Moreover, the 
Fig. 2 and Fig. 3 demonstrate the fact that the voids are mostly found in the 
neighborhood of the centers of anisotropy. Further, to strengthen our claim 
that the shape anisotropy induces voids, we compared two membranes. In one 
case we added rod impurities made up of two disks (Rod2) in the assembly 
of circular disks, and in the other case we added circular impurities of larger 
size, which occupied the same area as that of Rod2. We found that the former, 
being more anisotropic, induced larger and more numerous voids as compared 
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to the later, though they occupied the same area. 

Thus, reduced to the bare essentials, the anisotropy in the interaction potential 
of the constituents, is seen to be responsible for the induction of large voids. If 
studied from the perspective of energy minimization, as the potential becomes 
direction dependent, some positions of the constituents are preferred over the 
other positions. This induces large voids. 

These features show a remarkable similarity with the observations reported 
in certain biological experiments [16]. These experiments deal with the size- 
dependent permeation of non-electrolytes across biological membranes. The 
effect of doping on the permeation of large molecules was studied in these 
experiments. The liposome-membrane used in these experiments was made 
up of mixture of two types of lipids (cardiolipins and phosphatidylcholine) 
in a proportion 1:10. The understanding of the enhancement of transport in 
doped membranes needed an algorithmic statement. The ingredients at the 
algorithmic level involved: 

(1) consideration of the structure as a strictly 2-dimensional assembly 

(2) the cross sections of molecules being considered as constituents 

(3) interactions of the constituents via the Lennard Jones potential 

(4) permeating particles being considered as hard disks. 

The features reported in [16] bear a similarity with the simulation carried 
out with Rod2 as dopants. We have already seen in numerical simulations 
(See Fig. 5, curves (a) and (b)) that the Rod2 type of impurities induced large 
voids in the membrane. The appearance of larger voids naturally enhances the 
transport of large particles. Thus an enhancement in the transport of large 
non-electrolytes like glucose, which was observed in the lipid mixture [16] can 
possibly be understood using our simple approach. 

Further, apart from the biological implications, the model discussed is gen- 
eral enough to incorporate the studies of transport in various weakly bound 
granular media. 



5 Summary 

We have presented a numerical algorithm to compute the entire void statistics 
in a two dimensional membrane consisting of circular disks and dopants. We 
found that our simple two dimensional model has shown results consistent with 
features observed in a complex biological system. The biological justification 
of the model and implications are discussed elsewhere [17]. Nevertheless, our 
model and the proposed numerical algorithm which finds out the void statistics 
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in the model system are quite general and use no specific features of any 
particular system. Therefore it is possible to use this method effectively in 
various systems from diverse disciplines. The result that the shape anisotropy 
induces large voids in mixtures may be used as a tool for achieving controlled 
selective permeability across such a system by merely changing the shape of 
the constituents of the mixture. 

Acknowledgments We thank N.V. Joshi, Deepak Dhar, H. E. Stanley and 
S.S. Manna for fruitful discussions. 
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Figure captions: 

Fig. 1 Some examples of the impurities. 

(a) Rod type impurity made up of three circles (Rod 3 ). 

(b) L type impurity made up of four circles (L 4 ). 

(c) X type impurity made up of five circles (X 5 ). 

(d) T type impurity made up of five circles (T5). 

Fig. 2 Typical equilibrium configuration of a membrane without impurity. 

There are 556 circular disks used to form this membrane. The number is so 

chosen that the area occupied is ~ 70%. The a in Lennard- Jones potential 

is chosen as two times the radius of a circular disk. 
Fig. 3 Typical equilibrium configuration of a membrane with impurity of type 

Rods. The amount of impurity is 1:10 proportion. All the other parameters 

are same as Fig. 2. 

Fig. 4 Figure describes the algorithm to compute void statistics. The radius 

of a circular disks (black disks) is R. These disks are expanded by amount 

r, so that the region V is the void for particle of size r. 
Fig. 5 The graphs of number of r-voids as a function of r measured in units 

of the radius of the constituents. 

Curve a The void distribution without impurities. 

Curve b The void distribution with impurity of type Rod 2 . 

Curve c The void distribution with impurity of type L 4 . 

Curve d The void distribution with impurity of type Rodzi. 

Typically 10000 Monte Carlo steps are thrown away as thermalisation, and 

it is ensured that total energy is minimized. The curves are averaged over 

100 Monte Carlo steps. 
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Fig. 3 
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